Competition and cooperation in one-dimensional stepping stone models 
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Cooperative mutualism is a major force driving evolution and sustaining ecosystems. Although 
the importance of spatial degrees of freedom and number fluctuations is well-known, their effects 
on mutualism are not fully understood. With range expansions of microbes in mind, we show that, 
even when mutualism confers a distinct selective advantage, it persists only in populations with high 
density and frequent migrations. When these parameters are reduced, mutualism is generically lost 
via a directed percolation process, with a phase diagram strongly influenced by an exceptional DP2 
transition. 
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Cooperation is at the heart of many complex sys- 
tems [1, 2). On an organism level, gut bacteria help their 
hosts digest cellulose. On an ecosystem level, plants often 
rely on fungi to receive important nutrients. Even human 
societies are products of cooperation between individu- 
als. Despite the apparent advantage and pervasiveness of 
mutualistic interactions, their existence is often difficult 
to explain by a naive application of Darwinian natural 
selection: Cooperation can succumb to cheating [1] and, 
as we show here, to number fluctuations. 

To model complex interactions between individuals or 
species, Maynard Smith developed evolutionary game 
theory [2]. The central idea of game theory is that the 
fitness of an organism depends on the frequency and 
types of encounters with other organisms in the popu- 
lation. Evolutionary games are usually analyzed using 
mean- field- type approximations, which neglect both spa- 
tial correlations and number fluctuations. However, these 
simplifications are not appropriate for natural popula- 
tions living in spatially extended habitats and can miss 
important stochastic aspects of population dynamics. In 
particular, the interplay of stochasticity and spatial de- 
grees of freedom leads to spatial demixing of different 
species or genotypes in the population [3, 4], which can 
significantly decrease the probability of mutualistic inter- 
actions. 

Following pioneering work of Nowak and May [5] , sev- 
eral studies have investigated the effects of space on evo- 
lutionary games [6, 7] (and references therein) using sim- 
ulations on a two-dimensional lattice with a single non- 
motile individual per site. Although these studies under- 
scored the significant effects of spatial structure on evolu- 
tionary dynamics, outstanding issues remain. First, the 
outcomes of these lattice simulations are very sensitive 
to the exact rules of birth and death updates and inter- 
action pattern between nearest neighbors [7]; as a result, 
these studies do not smoothly connect with the well- 
understood dynamics in spatially homogeneous (well- 



mixed) populations. Moreover, it is not clear whether 
a model with a single nonmotile organism per site and 
nearest neighbor interactions is a good description of any 
species. Second, such models do not allow systematic in- 
vestigation of the role of migration and the magnitude of 
number fluctuations, which are important for the appli- 
cations of the theory to natural and experimental popu- 
lations. Third, closely related voter models in two spa- 
tial dimensions have very slow logarithmic coarsening [8] ; 
2d simulations typically do not explore the time scales on 
which spatial demixing of species becomes important. 

This letter studies a stepping stone model of one- 
dimensional population genetics [9]. The stepping stone 
model preserves the interaction pattern of well-mixed 
populations, but includes migrations as well as num- 
ber fluctuations, which are controlled by the population 
density. We focus on a one-dimensional model because 
stochastic effects are more pronounced in lower spatial 
dimensions [4]. More importantly, the spread of mu- 
tualism in two dimensions often occurs via a traveling 
reaction-diffusion wave, where the most important dy- 
namics often occurs at a moving quasi-one-dimensional 
frontier [4]. We find that, for one-dimensional popula- 
tions, mutualism persists in a much smaller region of 
parameter space than for well- mixed populations. Mu- 
tualism is particularly unstable against spatial demixing 
when the benefits to the interactants are unequal. The 
critical strength of mutualism required to sustain cooper- 
ation increases with migration rate and population den- 
sity. As the strength of mutualism is reduced, the pop- 
ulation undergoes a nonequilibrium phase transition in 
the universality class of either directed percolation (DP) 
or Z2 symmetric directed percolation (DP2); see Ref. [10] 
for a comprehensive review of DP models. 

The stepping stone model [9] we use to simulate 
the evolution of mutualism on a computer consists of 
demes (islands) arranged on a line, with spacing a. Each 
deme has N organisms, which can reproduce and mi- 
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grate. One reproduction and one migration update in 
each deme constitute a time step. One generation, in 
this Moran process [4], corresponds to TV time steps be- 
cause every individual in a deme is updated once on av- 
erage. Organisms can migrate to one of the two nearest 
neighbors with equal probability, and each organism has 
a probability m to migrate during a generation time r. 
Reproduction occurs within a deme by selecting a ran- 
dom individual to die and another individual to repro- 
duce. The probability to be selected for reproduction is 
proportional to individual's fitness. To study mutualism, 
we assume that the fitness is a sum of two contributions: 
a background reproduction rate, which we scale to one 
for all organisms in the population, and a benefit due to 
mutualistic interactions with other organisms in a deme 
(e.g. due to exchanging nutrients). Let the benefit to the 
organism of type i from interacting with the organism of 
type j be a^. If the types fractions within a deme are 
then the corresponding fitnesses W{ in a given generation 
are W{ = 1 + ctijfj because the increases in growth 
rate due to mutualism should be weighted by the density 
of cooperating organisms. 

For simplicity, we consider only two cooperating 
species (or genotypes) and let the frequency of species 1 
be x), where t is time, and x is position. The fre- 
quency of the other species is then 1 — f(t,x). In the 
limit of weak selection, when C 1, we find a contin- 
uum description of this one-dimensional stepping stone 
model in terms of a generalized stochastic Fisher equa- 
tion [4] 



% = sf(l -/)(/*-/) + D s g + y/D g f(l-f)r(t, x) , 

(1) 

where T(t,x) is an ltd delta-correlated Gaussian white 
noise 1 , D s is the spatial diffusion constant, and D g ~ 
1/N is the strength of number fluctuations. The key 
parameters s = (ai + a 2 )/r and /* = a±/s are given 
by ai = (ai2 - a 22 )/r and a 2 = (a 21 - on)/r. The 
selective advantage (or strength) of mutualism is given 
by 5, while /* is the equilibrium fraction of species 1 
that would occur in a spatially homogeneous population 
without number fluctuations. 

The usual mean-field treatment neglects spatial corre- 
lations and fluctuations. With the neglect of the last two 
terms, Eq. (1) becomes an ordinary differential equation, 
and its dynamics can be easily analyzed. There are four 
possible outcomes as shown in Fig. 3a and the Appendix. 
The population develops mutualism when a\ and a 2 > 0, 
one of the species outcompetes the other when aia 2 < 0, 
and, when ol\ and a 2 < 0, the population is bistable, with 



1 In Ito's formulation, f(to, x) and T(to, x) are independent for any 
given to, but special rules of Ito calculus must be used to differ- 
entiate a composite function; see Ref. [4] for further discussion. 



the either species capable of outcompeting the other de- 
pending on the initial conditions. 

We first explore number fluctuations in a population 
without spatial structure, e.g. consisting of one deme or 
a finite population of 7V tot organisms with very large D s . 
A stochastic treatment must also account for absorbing 
boundary conditions at / = and / = 1, when one of the 
two species goes extinct. The absorbing boundaries arise 
because there is a finite probability to find the population 
in any of its discrete states when the population size is 
finite. Therefore, after a sufficiently long time, the popu- 
lation will reach one of the absorbing boundaries and be- 
come fixed. The splitting probabilities and fixation times 
can be calculated for this zero-dimensional problem with 
fluctuations using the Kolmogorov backward equations; 
see the Appendix. 

In a spatially extended population, however, local ex- 
tinctions can be prevented or rescued through migration. 
Suppose, in particular, that migrations are frequent and 
mutualism is sufficiently strong to keep the population 
near the equilibrium fraction /*. In this limit, we can 
extend the mean-field approximation to account for fluc- 
tuations and spatial degrees of freedom by replacing the 
nonlinear reaction term in Eq. (1) with a linear one: 



dt 
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f)+D S 0± + y/D g f(l-f)T(t,x). 

If f(t,x) ~ /*, the error we make should be small; more 
importantly, Eq. (2) can now be solved exactly. The solu- 
tion is most easily obtained in terms of the average spatial 
heterozygosity H(t,x), a two-point correlation function 
equal to the probability to sample two different species 
distance x apart: 

H(t, x) = {f(t, 0) [1 - f{t, x)} + fit, x) [1 - fit, 0)]) . (3) 

Using the Ito calculus, we derive the equation of motion 
for H{t,x) from Eq. (2), 



OH 
~dt 



2D S |^ - D g 5{x) - sH*(l - H*) 



H, (4) 



where S(x) is the delta function, H* = 2/*(l - /*), and 
we, for simplicity, assumed that the population is initially 
uniform with the equilibrium fraction /(0, x) = /*. The 
stationary solution, valid at long times, reads 



i^(oo, x) 
H* 



= 1 



V^8f*(l-f*)/D a 



l + j&8D a f*(l-f*)/D* 



(5) 



Since H* is the heterozygosity of a well-mixed population 
with / = /*, the fraction on the right hand side is the 
correction to the mean- field analysis. Thus, we see that, 
for s <C D^/D S1 the probability of the two species coex- 
isting at any particular point in space [given by iJ(oo, 0)] 
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FIG. 1: (Color online) Mutualism in the one-dimensional stepping stone model, (a) Spatial demixing for N — 30, mN — 1, and 
no interspecies interactions, all a%j — 0. Green (light gray) and red (dark gray) represent species 1 and 2 respectively. Every 
deme and every tenth generation are shown, (b) The same as in a, but with strong mutualism an — ol2\ — 0.5. (c) Heat map 
of H(4 • 10 6 , 0) from simulations with the same parameters as in a, but with 10 4 demes and varying a\2 and 021. 



becomes small, which is inconsistent with mutualism and 
our assumption that /(£, x) ~ /*. Hence, we anticipate a 
critical value of s below which mutualism must give way 
to spatial demixing. 

Although the hierarchy of moment equations does not 
close for the original, nonlinear problem given by Eq. (1), 
the average spatial heterozygosity H(t,x) is still useful 
for characterizing the behavior of the system. In particu- 
lar, the average local heterozygosity H (£, 0) can be used 
to measure the amount of mutualism. Equation (5) sug- 
gests that i7(£, 0) reaches a nonzero steady state value 
when s ^> D^/D s . However, when s = 0, the exact solu- 
tion of Eq. (1) reveals that instead of reaching a steady 
state, H(t, 0) decays to zero as t~ x l 2 [4]. When species do 
not coexist locally, mutualism is impossible. Hence, we 
can use the long time behavior of H (£, 0) to distinguish 
between populations where mutualism can and cannot 
persist. See the Appendix for another quantity to distin- 
guish the phases, similar to the susceptibility in equilib- 
rium physics. 

The phase digram obtained from simulations is shown 
in Fig. lc. The region of parameters where mutualism 
can evolve is significantly reduced compared to the well- 
mixed prediction shown in Fig 3. In particular, mutual- 
ism is impossible even for positive s, provided s is small. 
Fluctuations and spatial structure also favor symmetric 
mutualism, with a\ « oli, i-e. when the two species 
benefit equally from the interaction. The mutualistic 
phase [characterized by lim^oo H(t, 0) 7^ 0] is separated 
from the demixed phase [lim^oo H(t, 0) = 0] by two lines 
of second order phase transitions that meet in a cusp: 
lim^oo H (£, 0) decreases continuously to zero as these 
lines are approached. 

Nonequilibrium phase transitions from an ac- 
tive (mixed) to an absorbing (one of the two species) 
state have been studies extensively; see Ref. [10]. 
Generically, when the absorbing states are not symmet- 
ric, ct\ ^ O62, the exit from mutualism belongs to the DP 
universality class. We can most readily see this for /* 



close to an absorbing boundary, say f « 1. For large 
mutualistic selective advantage s, species 1 then remains 
at low frequencies throughout the population. As the 
strength of mutualism is decreased, some spatial regions 
stochastically lose species 1, but the more abundant 
species 2 persists. Local extinctions are opposed by 
the spread of species 1 from the nearby regions via 
Fisher waves. This dynamics is just that of DP in [10]. 
When a\ = the absorbing states are symmetric 
and the local extinctions of either species are equally 
likely. As a result, this phase transition belongs to DP2 
universality class. We checked that our simulations are 
consistent with the DP2 "bicritical point" by calcu- 
lating how H(t,0) decays for different values of s in a 
population that is initially well-mixed (see Fig. 2a) and 
then collapsing these decay curves onto a unique scaling 
function using DP2 exponents as shown in Fig. 2b. 
Equation (1) is also known to describe a DP2 transition 
for /* = 1/2 [11]. Although the DP2 transition occurs 
only at a point, it influences a large portion of the phase 
diagram and governs the nonlinear shape of the DP 
transition lines near this "bicritical point." 

To understand how phase boundaries depend on the 
parameters of the model, it is convenient to measure 
distance in the units of D s /D g and time in the units 
of D s /Dg. When, Eq. (1) is nondimensionalized, and 
the dynamics is controlled by only two dimensionless pa- 
rameters, /* and s = sD s /D^. We confirm this data 
collapse in simulations, see Fig. 2c. 

Spatial structure and number fluctuations change not 
only the mutualistic region (ai,a2 > 0), but also the 
whole phase diagram for well- mixed populations, see 
Fig 3. In particular, there is no bistable phase (with 
splitting probabilities sensitive to the initial conditions) 
in Id spatial populations. For almost all initial condi- 
tions, domains of species 1 or species 2 appear because 
of number fluctuations; the subsequent behavior can be 
analyzed in terms of the Fisher wave velocities of the do- 
main boundaries. For ol\ > ct2, this velocity is directed 
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FIG. 2: (Color online) Properties of the DP2 phase transition (a± — 0^2) from simulations with the same parameters as in 
Fig. lc. (a) Decay of the local heterozygosity with time; lower values of the mutualism selective advantage s lead to faster 
decay, (b) The collapse of if (£, 0) onto a universal scaling function shown in the inset with DP2 critical exponents for different s 
above the phase transition at s c — 0.109. For the DP2 critical point, we expect /3 = 0.92 and v\\ — 3.22 [10]. (c) H{t — 4- 10 6 , 0) 
as a function of s — sD s /D 2 g for different values of m and N. We show how to obtain D s /D 2 g in the Appendix. 



from species 1 to species 2, and the direction is reversed 
for a\ < ct2. As a result, one of the species takes over, 
much like an equilibrium first order phase transition pro- 
ceeds through nucleation and growth. 

When ol\ = 0x2 and mutualism is unstable, population 
segregates into single species domains, and the dynamics 
is driven by the random walks of domain boundaries. 
For s > this demixing is slowed down by mutualism, 
but for s < it is sped up initially due to the reaction 
term in Eq. (1). After domains form, however, more 
negative values of s lead to slower coarsening because 
the diffusion constant of domain boundaries decreases. 
Surprisingly, the exactly solvable limit of a\ = = 
undergoes the fastest demixing in the long time limit for 
large system sizes (see Fig. 4). 

We have shown that a critical strength of mutualism 
is required to overcome the spatial demixing of species 
driven by local number fluctuations. The critical strength 
of mutualism strongly depends on the symmetry of the 
interaction, and mutualism is more likely to evolve be- 
tween species that share the benefits equally. We believe 
that these predictions can be tested in growing microbial 
colonies on a Petri dish because colony frontiers behave as 
quasi-one-dimensional populations [4] and because mutu- 
alistic interactions can be engineered in the lab [12]. Such 
experimental studies would be of value not only for evo- 
lutionary biology, but also for nonequilibrium statistical 
mechanics because they could provide important experi- 
mental tests of the theory. 

After submitting this paper, we learned of a preprint 
by Dair Astra et a/., "Strong noise effects in one- 
dimensional neutral populations" (arXiv: 1012. 1209), 
which discusses symmetric cooperation in a similar 
model, corresponding to /* = 1/2 and a\ = OL2 in our ter- 
minology. Our work was constructed with experiments 
at microbial frontiers in mind; hence, it differs due to its 



focus on asymmetric interactions ol\ and c\2 of arbitrary 
sign and the large deme sizes in our simulations. 
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Appendix A: Supplementary Material 

Competition and cooperation in a zero-dimensional 
model 

Here, we present a more detailed analysis of evolution- 
ary game theory in a spatially homogeneous (well-mixed) 
population with N to t organisms. We consider only two 
species, so /(£), the frequency of species 1, is sufficient 
to fully characterize the state of the population. The 
equation of motion for f(t) takes the following form 

j t f{t) = 8f(t)[l-f(t)]\T-f(t)] + y/V g f(t)[l ~ f(t)]T(t), 

(Al) 

where T(t) is a zero mean Gaussian white noise inter- 
preted according to the ltd prescription. 

<r(ti)r(t 2 )> = (j(ti-t 2 ), (A2) 
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FIG. 3: (Color online) Phase diagrams of cooperation and competition, (a) The phase diagram in the mean-field limit, which 
neglects spatial structure and number fluctuations, (b) Schematic phase diagram for evolutionary games in one dimension, (c) 
Schematic plot of the change of species equilibrium frequencies along a cut a\ +«2 = const intersecting the DP phase transition 
lines through the yellow region in the phase diagram shown in (b). The two dashed lines show the locations of the DP phase 
transitions. We find that the nonlinear dependence of the /i - /2 on ai - «2 near DP transitions is consistent with the DP 
exponent f3 = 0.276486 [10]. (d) Illustration of the fact that the phase behavior summarized in (a)-(c) is not specific to the 
well-mixed initial conditions used for most of our simulations. Green (light gray) and red (dark gray) represent species 1 and 2 
respectively. The simulation parameters correspond to the yellow, mutualistic phase in (b). In fact, the same parameters are 
used as in Fig. lb, but with the different initial conditions. At t = 0, the population consists of two equal domains: one of 
species 1 and one of species 2. Due to the periodic boundary conditions, this initial condition leads to two domain boundaries. 
The system rapidly forms a mutualistic phase, which squeezes out the red and green domains. 
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and V g ~ 2/(7V tot r) measures the strength of number 
fluctuations. 

In the absence of noise, the system has four generic 
types of behavior. 

• a\ > 0, 0L2 > corresponds to the snowdrift game 
or mutualism [1]. The stable fixed point is at / = 
/*, which is between zero and one. 

• ol\ < 0, d2 < corresponds to the coordination 
game. There are two stable fixed points / = 
and / = 1. /* is still between zero and one, but 
this fixed point is unstable. The long time fate of 
the system is entirely determined by whether the 
initial fraction of species 1 is above or below /*. 

• ql\ > 0, ol<2 < corresponds to competitive ex- 
clusion. Species 1 dominates, and the stable fixed 
point is at / = 1. 

• ol\ < 0, ot2 > also corresponds to competitive 
exclusion. Species 2 dominates, and the stable fixed 
point is at / = 0. 

In the presence of noise, the system always reaches 
one of the absorbing states at / = or / = 1 at long 
times. The splitting probabilities can be calculated from 
the backward Kolmogorov equation [13] , and we find 
that J7(/o), the probability to eventually reach / = 1 
starting from f = fo, is given by 



where P(t, /) is the probability distribution of / at time t. 
The Fokker-Plank equation is more intuitive for most 
physicists when the diffusion coefficient does not depend 
on / because, in this case, it describes diffusion of a par- 
ticle in a potential. We can achieve such a simplification 
by the following change of variables 



/ = sin 2 (p/2). 
The potential is then given by 



(A6) 



V = V g In[sin(p)] + - cos(p) [f* ~ -J + - cos 2 (p). (A7) 

For /* = 1/2, V has a minimum in the interior only 
for s > AV g suggesting that the decay of weak mutualism 
is not protected by a barrier, which is also true for other 
values of /*. 



Connection between H(t,x) and a susceptibility 

In the main text, we use lim^oo H (t, 0) to distin- 
guish between mutualistic and non-mutualistic phases. 
Although this quantity is very convenient and easy to 
measure in simulations, it is rarely used in equilibrium 
physics to characterize phase transitions. We can make 
a more explicit connection between the equilibrium the- 
ory of phase transitions and the out of equilibrium loss 
of mutualism by considering the following quantity 



U(fo) 



f 3 
Jo 



df 



He*"-™* 



(A3) 



In the limit of strong noise V g ^> s, the fixation proba- 
bility is approximately given by 



U(fo) = Ml + 3^(1 - 3/* + 3/ - /o 2 )]- (A4) 

Note that for / <C /* we recover the 1/3-law [14], which 
states that, for s < the fixation probability of species 1 
is greater than it would be under purely stochastic (neu- 
tral) dynamics, provided /* > 1/3. 

The average time to reach fixation can be calculated; 
see [13]. Cremer, Reichenbach, and Frey [15] have re- 
cently carried out the analysis of fixation time and found 
that it scales as 1/V g for \s\ <C V gi but becomes expo- 
nentially large when mutualism is strong. 

Another way to analyze the effects of noise on the dy- 
namics is to look at the Fokker-Plank equation corre- 
sponding to Eq. (Al), namely 



_9 

&t 



p(tj) = -s^[f(i-f)(r-f)p(tj)} + 



(A5) 



x(t) 



f 

J — c 



[H(t,x) - lim H(t,x)]dx. 



(A8) 



Note that H(t,x) < linx^oo H(t, x), so x 18 positive 
and measures the degree of species demixing. Since x 
is an integral of a two-point correlation function H(t,x), 
it is analogous to the susceptibility used in equilibrium 
physics. 

In the mutualistic phase, both species are present, 
so lim t ^oo lim^oo H(t, x) > 0. Then lim^oo X (*) is fi- 
nite and measures the degree and spatial extent of the 
reduction of the species diversity due to number fluc- 
tuations. Close to one of the phase transition lines in 
Fig. 3b, lim^oo H (t, 0) « 0; therefore, lim^^ xif) ~ 
£ lim^^oo H(t, x), where £ is the correlation length. As 
the phase transition is approached, £ diverges with an ex- 
ponent z/j_, which is 1.096854 and 1.83 for DP and DP2 
transition respectively. In addition, lim^oo H (t, x) = 
1/2 for the DP2 transition, but linx^oo H(t,x) vanishes 
with exponent f3 = 0.276486 for the DP transition [10]. 

In a non-mutualistic phase, x(t) diverges with time be- 
cause H (t, x) « for x smaller than the typical domain 
size, which grows in time due to coarsening. Close to the 
phase transition this divergence occurs with the expo- 
nent 1/z, where the dynamical exponent z = 1.580745 
and z = 1.74 for DP and DP2 transitions respec- 
tively [10]. For the exactly solvable case of s = 0, x(t) ~ 
t 1 ' 2 . 
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Estimating D g and D s 

Generalized stochastic Fisher-Kolmogorov equations, 
e.g. Eq. (1), describe population dynamics on long time 
and length scales; therefore, the diffusion constant D s 
and the strength of number fluctuations D g should be 
viewed as phenomenological parameters in an effective 
field theory. As a result, they may often be complicated 
functions of the microscopic parameters. For the step- 
ping stone model, D s = ma 2 /(2r) and D g = 2a/ (Nr) 
in the limit of large N and small m [4]. Our computer 
simulations, however, are not in this limit, and we have 
to estimate D s and D g numerically. 

This estimation can be easily done from simulations 
with all dij = because the stepping stone model without 
species interactions can be solved exactly. In particular, 
the decay of H(t,0) from a homogeneous 1:1 mixture of 
the species satisfies the following relation [4] 



H(t,0) 



2D, 



-1/2 



+ 0{r z ' 2 ). 



(A9) 



We can then obtain D s /D 2 g by fitting Eq. (A9) to the 
simulation results. One can also estimate D s and D g 
individually from the average size of domains; see Ref. [4]. 

In our letter, we used the aforementioned estimate 
of D s /D g to collapse H(t = oo,x = 0) as a function 
of s in Fig. 2c. We found that the DP2 transition occurs 
at s = 0.80. 



FIG. 4: (Color online) Decay of H(t, 0) in initially well-mixed 
populations with N — 30, mN — 1, an — a<2i — 0, a\i — a,2i, 
and 10 4 denies. Note that for more negative a±2 the initial 
decay is faster, but it crosses over to a slower decay at longer 
times. 
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